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Various problems in materials modeling can only be ad- 
dressed by "empirical" interatomic potentials^. Say we wish 
to evaluate a physical property (e.g. total energy) of some ma- 
terial with a complex crystal structure. We cannot directly in- 
sert the results of crystallographic refinements, as they (almost 
always) include sites with mixed or fractional occupancies. To 
obtain valid results, we must assign those occupancies plau- 
sibly based on computed energy differences. A single relax- 
ation with fast ab-initio codes such as VASP 2 may be tractable 
in a cell of 10 3 atoms; however, if (say) 5% of those atoms are 
uncertain, this must be repeated many times to find the one 
optimum state out of 2 50 possibilities. 

But with empirical, approximate potentials that can be eval- 
uated in negligible time, this optimization is tractable, and can 
be followed up by ab-initio relaxation to obtain the most ac- 
curate positions and total energies. When combined with ef- 
fective search algorithms, such as genetic algorithms or the 
"cell-constrained melt-quench" method (presented in Sec.|II| 
this is a powerful tool for ab-initio structure discovery. Some 
other questions that call for empirical potentials are phonon 
spectra (or other dynamics), and thermodynamic simulations 
of phase transformations in complex alloys. 

A popular framework of empirical potentials is the 
"embedded-atom method" (EAMJ2I4J, in which the full Hamil- 
tonian contains the usual pair term Vij(ri — tj), but also an 
implicitly many-atom term J2i U(p(ri)), where p(rj) is a sum 
of contributions at atom i from nearby atoms. Accurate EAM 
potentials are straightforward to extract for pure elements,but 
demand patience and skill to obtain even for binary systems; 
there is no systematic recipe for multicomponent systems. 

Here we present an alternative approach fitting only pair in- 
teractions but incorporating Friedel oscillations. These "em- 
pirical oscillating pair potentials" (EOPP) have the form 6 

V(r) = % + -^-cos(fc*r + &) (1) 

All six parameters, including fc*, are taken as independent in 
the fit for each pair of elements. Eq. ([T]i was inspired by ef- 
fective potentials (e.g. Refs.[7]and[8]) used in previous work 
on structurally complex metals, e.g. quasicrystals 9 10 . In such 
systems, energy differences between competing structures are 
often controlled by second- and third-neighbor wells due to 



Friedel oscillations, which are a consequence (mathemati- 
cally) of Fourier transforming the Fermi surface, or (physi- 
cally) are equivalent to the Hume-Rothery stabilization by en- 
hancing the strength of structure factors that hybridize states 
across the Fermi surface!^ In Pettifor's framework (Ref. [TT] 
Sec. 6.6), the short-range repulsion is captured by the first 
term of Q; the medium-range potential (first-neighbor well) 
as well as the long-range oscillatory tail are captured by the 
second term, their relative weights being adjusted by the r\\ 
and r\2 parameters. Since the second term has this double duty 
of representing both the nearest-neighbor and long-distance 
behaviors, the fitted 1 / r m decay generally does not agree with 
analytic asymptotic result of 1/r 3 ; also, fc* need not match the 
Fermi wavevector, and these parameters take different values 
for the six kinds of pair potentials. 

In the rest of this Letter, we describe how our potentials are 
fitted (typically for a particular composition range) and then 
demonstrate their capabilities through case studies in the alloy 
system Al-Cu-Sc, which has local order similar to the binary 
quasicrystal i(CaCd) or to Zn-rich Mg-Zn alloys. Along with 
this, we also describe a method of "constrained cell quench- 
ing" that accesses low-energy structures in surprisingly large 
cells. Finally, we will summarize other systems where EOPP 
have been applied, and discuss their limitations. 



I. DATABASE AND POTENTIAL FITTING 

The parameters in Eq. (fTl are fitted to an ab-initio dataset 
(using the VASP code ) combining both relaxed T = struc- 
tures and molecular dynamics (MD) simulations at high T. 
Four criteria in choosing structures for the database are (1) 
to bracket the composition range of interest; (2) to mix sim- 
ple and complex structures; (3) to ensure adequately many 
contacts of every kind (in particular, nearest neighbors of the 
least abundant species); (4) to have similar atom densities^. 
In ab-initio MD simulations of the simpler structures, a super- 
cell is always used with dimensions comparable to the fitting 
potential cutoff radius, which (for the fit procedure) is always 
12A . 

We define each structure's energy as a difference relative 
to a coexisting mixture (with the same total composition) of 
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FIG. 1: (a). Top: EOPP form (Eq. [Q fitted to "GPT" potentials 
(Ref.ll6t; the fitted/GPT curves lie exactly on the top of each other. 
The Cu-Sc potential is not available from Ref. 8; the figure shows a 
fit to ab-initio data with the other five pairs constrained to GPT form. 
Bottom: the EOPP form fitted to VASP force and energy data. The 
curves fitted to "GPT" are shifted by 0.15 eV/atom for clarity, (b). 
Scatter plot of the forces for the EOPP fit to [Eq. {T}], from the fitted 
Al-Cu-Sc potentials shown in (a), bottom. The pair potential forces 
Fj are vertical axis, ab-initio data horizontal. 



Our empirical potentials can be compared with Moriarty's 
"GPT" potentials 16 , derived from a systematic expansion, 
which are known for all but one of the six Al-Cu-Sc pair 
potentials [see Fig. |TJb)] . First, when the GPT potentials are 
fitted to the EOPP form (Eq. ([T]l, the curves virtually lie on top 
of each other: thus, the EOPP form can represent all features 
of the GPT potentials. Second, the GPT and EOPP potentials 
show a strong similarity; the main mismatches are (i) EOPP 
has no first-neighbor well for Al-Cu and Cu-Cu (ii) the os- 
cillation wavelength for Sc-Sc differs (iii) the overall energy 
scale of the GPT pair potential is too small by ^50% In 
practice empirically fitted potentials account for some of the 
many-body contributions by modifying their pair terms, and 
hence work better than truncating a systematic expansion (e.g. 
GPT 16 ) after the pair terms. 



II. CONSTRAINED-CELL MELT QUENCHING 



reference phases, chosen to bracket all database compositions. 
Every structure is used for both forces (from MD at high T) 
and energy differences 13 (high-T MD, as well as relaxed at 
T = 0). For the high-T portion, we took one snapshot of 
each structure at the end of a short ab-initio MD run. Usually 
> 10 3 force components enter the fit, along with ~ 50 energy 
differences. Typical forces are ~ 3 eV/A and typical energy 
differences are ~ 0.3 eV/atom; the fit residuals are ~ 5% and 
~ 1% respectively [see e.g. Fig.[TJa)]. 

Our least-squares fit minimizes (by the Levenberg- 
Marquardt algorithm) x 2 = £ AEf/a E 2 + £ |AF,-| 2 /cr F 2 , 
where {AEi} and {AF 3 } are the energy and force residuals; 
we found a weighting ratio oeI^f ~ 10 _3 A was optimal so 
that neither energies nor forces dominate the fit. 

There is some risk of converging to a false minimum (or not 
converging at all, from an unreasonable initial guess). Thus, 
it is important to repeat the fit from several starting guesses. 
For this we used, e.g., potentials first fitted to pure elements 
or binary systems, and also used a library of parameter sets 
previously fitted for some different alloy system. The fitted 
parameters in ([TJ for each of our examples are available as 
supplementary information 6 ; similar potentials were plotted 
in Refs.[H(for Al-Mg) and [BJ (for Sc-Zn). 

For our specific case of Al-Cu-Sc, our database had 84 
relaxed energies at T = OK from (besides the pure ele- 
ments) the binaries Al 3 Sc.cP4, Cu 2 Sc.tI6, Al 2 Sc.cF24, and 
Al 2 Cu.cF12 and ternaries AlCuSc.hP12, Mg 2 Cu 6 Ga 5 .cP39 
(with Mg-s-Sc, Ga— >A1), and AlCrCu 2 .cF16 (with Cr-^Sc); 
structures are identified by Pearson labels. The database 
also had 7428 force points at T > OK, taken from all 
those structures, and additionally from Al 2 Sc.hP12, and 
Sc(Al 1 _ 2 .Cu :I .)g.cI168. The EOPP potentials give very good 
agreement with this database, as shown in Figure [TJa): the 
r.m.s. deviation of forces was 0.11 eV/A and that of energy 
differences was 21.7 meV/atom, with relative weights set to 
maximize accuracy of the force data. 5 Values of the parame- 
ters are listed in Table|I] Appendix [A] 



We now turn to the second method which, together with 
fitted potentials, has enabled the present structure study and 
others. In many alloy systems, with no structural informa- 
tion known except the composition and the unit cell, struc- 
tures with > 100 atoms per cell may be predicted from care- 
ful "melt quenching" (MQ) simulations. The relation to the 
EOPP notion is that the larger systems - particularly when 
supercells are used - are too large for direct ab-initio calcula- 
tions, so empirical potentials are crucial for melt-quenching. 
This method has been applied with GPT potentials to improve 
known structures of the decagonal Al-Co-Ni quasicrystals 18 . 

In most cases, annealing requires "tempering Monte Carlo 
(MC)'G21 wherein ~ 10 samples are annealed simultaneously 
at equally spaced temperatures spanning across the melting 
temperature. Each tempering cycle includes a short MD run 
(~ 1 ps with 1 fs steps) followed by lattice-gas MC annealing 
(~ 200 attempts per atom) in which the chemical identities of 
two randomly selected atoms may be swapped. Then, pairs 
of samples may be swapped using a Metropolis-like criterion. 
(For the Al-Cu-Sc structures we studied, our temperatures 
spanned 600-1700 K with spacing AT =100K, and the inter- 
action cutoff was set to r cut = 9A.) The resulting structure 
may be tested subsequently by ab-initio calculations of the to- 
tal energy, for which we used VASfl 

A key diagnostic in a tempering simulation is the (time- 
dependent) energies E m (t) of all samples, as specific heat is 
approximately AE/AT, where AE — E„ l+ i — E m . For 
our structures AE peaked around T =1400K, indicating the 
melting point. In our cases with 52 or 84 atoms per primitive 
cell, several independently initialized runs yielded identical 
low-T structures, which took 100-200 cycles. 



III. TARGET STRUCTURES 

We now apply empirical pair potentials and cell- 
constrained melt-quenching to realistic structure predicti on in 
the Al-Cu-Sc system, for two newly discovered phases^El 
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in which lattice parameters were known from electron mi- 
croscopy but no single-grain structures were available for 
structure determination. We label each phase by its Pear- 
son symbol. These belong to a family of phases in which 
atom size is the salient attribute: Ca, Mg, Sc, and rare earths 
are "large" atoms, constituting ~ 1/6 by number; whereas Al, 
Zn, Cu, and Cd are the majority "small" atoms. 

Our first target is "c/168" with 84 atoms per primitive cell, 
hypothesized to be isostructural with the ScZng phase. No 
single-crystal data are available; a preliminary Rietveld refine- 
ment of powder data technique 21 confirmed the ScZng type 
structure, with refined lattice parameter a=13.52A. The chem- 
ical ordering and occupancies of Al/Cu could not be deter- 
mined reliably prior to our modeling 2 ^. 

The ScZn6 phase is an "approximant" of the recently dis- 
covered thermodynamically stable, icosahedral quasicrystal 
z(ScZn), meaning its unit cell contents are identical to a frag- 
ment of the quasicrystal structure. i(ScZn) has the same 
atomic structure as i(CaCd) 22 , which (along with similar rare 
earth-Cd quasicrystals) are the only known stable binary qua- 
sicrystals. Each site is specific to either a small atom (Cd or 
Zn) or a large atom (Sc, Ca, or rare earth). 

Both c/168-ScZn6 and the related quasicrystals are under- 
stood to be packings of a three-shell, icosahedrally symmet- 
ric cluster called the "Tsai cluster". Its outermost shell has 
an icosahedron of 12 large atoms plus 30 small atoms (form- 
ing an icosidodecahedron) on the midpoints of the large-large 
bonds [See Fig.[2|a)], and inside that is a dodecahedron of 20 
small atoms. At the cluster center is a tetrahedron of small 
atoms; as this breaks the icosahedral symmetry, it has many 
degenerate orientations, which leads to interesting slow dy- 
namics of reorientations^ and to structural ordering transi- 
tions at low T in the related crystals. 

Our second target is "oCT04", of composition 
AI38.8CU45.7SC15.5, Initially ED the lattice parameters 
a=8.32A, 6=8.36A, c=21.99A and C-centered orthorhombic 
Bravais lattice were identified from powder data but that was 
insufficient to refine the structure. 



IV. RESULTS: PHASE STABILITY AND PSEUDO-TSAI 
CLUSTERS 

Constrained-cell melt-quenching using the same database- 
fitted EOPP potentials converged to low-temperature struc- 
tures in both cases, the energies of which were subsequently 
evaluated by VASP. In the c/168 case this was the known 
structure. In the oC104 case, a previously unknown struc- 
ture was obtainecP, computed unstable by only 5 meV/atom 
with respect to c/168; its validity was confirmed by a subse- 
quent refinement of the powder data 20 with this as the starting 
point. This predicted structure had space group Amm% in 
which there were 1 1 Cu, 7 Al and 4 Sc Wyckoff sites, with no 
chemical disorder. 

The fundamental motif in both c/168 and oC104 is an 
icosahedral cluster whose innermost shell has less symme- 
try. In c/168 this is a Tsai cluster [Fig. [2] (a)] in which each 
of the three Al/Cu shells has composition AI0.5CU0.5. Strik- 




FIG. 2: (a) Tsai cluster in c/168 structure (left) and (b) "pseudo- 
Tsai" cluster (right) found in the oC104 structure. View along 2- 
fold direction. In each panel, the successive shells are (1) tetrahe- 
dron AI2CU2 or octahedron Cu6 (2) dodecahedron (Al,Cu)2o, and 
(3) icosidodecahedron (Al,Cu)3o plus icosahedron SC12. Large white 
circles represent Sc, dark smaller circles Cu, and light smaller circles 
Al. 



ingly, within each shell, the Al and Cu are segregated into 
hemispheres centered on fivefold axis of the icosahedrorl^, 
furthermore, the Al parts of each shell overlay the Cu parts 
of the preceding one and vice versa [see Figure |2|a)] . Due 
to this symmetry breaking, one expects the low-temperature 
phase to have at least several equivalent domains, correspond- 
ing to different orientations of the pseudo-Tsai clusters. 

In the oCT04 case we find a variant motiPS that we call 
the "pseudo-Tsai" cluster, in which the innermost shell is a 
Cu6 octahedron [Fig |2fb)] . This cluster was first described in 
Mg 2 Zn 11 .cP39 2 ^, where the outer shell is a Zn 8 Mg 12 dodec- 
ahedron plus a Zni2 icosahedrorPn 

The experimental (powder data) refinement 2 ^ of oC104 
is an average structure, with higher symmetry (space group 
Cmmm) than our model, accomodating 3 Sc Wyckoff sites, 
and 13 mixed Al/Cu sites (one of the latter being half- 
occupied). Each experimental site derives from one or two 
sites in our model, displaced on average by -0.17A (Al), 
~0.12A(Cu) and 0.08A(Sc). Crystallographic data of the 
low-temperature model structure, are compared with the ex- 
perimental structure in Appendix [B| The mean deviation of 
the refined Al/Cu content from the averaged model content is 
16%. 

The electronic density of states for c/168 has a deep, nar- 
row pseudogap at the Fermi energy, which tends to stabilize a 
unique composition. In contrast, oC104 has a shallower and 
wider pseudogap, suggesting a range of degenerate composi- 
tions, so that substitutional entropy might stabilize this phase 
at higher temperatures. 

The pseudo-Tsai clusters in oC104 (and also the less com- 
plex cP39) adjoin by sharing atoms such that their centers 
are closer by a factor ~ 1.618 (golden ratio) than Tsai clus- 
ters would be. Could pseudo-Tsai clusters be the basis of the 
newly reported i(AlCuSc) quasicrystal 30 ? 



V. CONCLUSION. 

We have shown that empirical potentials with the simple 
oscillating form ([T]), fitted from ab-initio data and combined 
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with a "cell constrained" brute-force quenching, allows de- 
tailed predictions of fairly complex low-temperature optimal 
structures in Al-Cu-Sc alloys, based on the very limited input 
of known lattice parameters and composition. Finding the cor- 
rect structure depends sensitively on having a quantitatively 
realistic potential, which is achievable only if that potential is 
constructed or fitted from ab-initio calculations. One would 
expect that the oscillating analytical form ([TJ is natural only 
for simple metals (e.g. Al or Mg, for which it does very well). 
But in fact, the EOP potentials sometimes work quite well 
even when angular or many-body interactions are important, 
e.g. the transition metal neighbors in Al-Cu-Sc. But - not 
surprisingly - they do poorly for elemental Zn or Ga. Of 
course, any pair potential fails when the electron density has 
large variations in space (as at vacancies, edge dislocations, 
or surfaces), or in the (many) cases where bond directionality 
(due to covalent bonding) is prominent. 

We believe the EOPP potentials are quite broadly applica- 
ble to mimic the atomic interactions of many metallic systems 
with sufficient accuracy to stand in for ab-initio energies when 
those would be computationally prohibitive. Although the 
EOP potentials were not formally presented before this paper, 
they have already been applied to a variety of intermetallics: 
(1) the site contents in complex structures, e.g. Al-Mg^or 
Al-Zn-Mg 32 ; (2) solving the complete structures of complex 
Mg-rich Mg-Pd phases (with > 400 atoms/cell) in conjunc- 
tion with diffraction, when the latter alone is insufficient 34 ; 
(3) phonon spectra, e.g. in Zn-S<JE and Mg-Zn 26 alloys, and 
even in liquid Bi-Li alloys 31 ; (4) the dynamics of the Tsai clus- 
ter tetrahedra in the c/168 structure ScZn 6 23 , as well as the ar- 
rangements of the asymmetric inner Alio shell in the pseudo- 
Mackay icosahedral clusters in quasicrystal-related Al-Ir, Al- 
Pd-Mn, and Al-Cu-Fe phased^]. 
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Appendix A: EOPP parameters for Al-Cu-Sc system. 

The EOPP potential parameters fitting Eq[T]to DFT/VASP 
data are shown in Table U 



Appendix B: Crystallographic data of the oC104 structure 

Table lists Wyckoff orbits of the low-temperature struc- 
ture, resulting from the melt-quenching procedure described 
in Sec. |n] using EOPP Al-Cu-Sc potentials (Tab. |IJ. The 
structure was subsequently optimized using DFT code VASP. 
The experimental (relaxed) cell parameters were a=8.3370 
(8.3247) A, 6=22.0150 (21.868) A and c=8.3370 (8.3247) A, 
space group Amm2 (no. 38). Note that a and c axis are 
swapped with respect to the setting in Ref. 20 Column (i 
is site multiplicity, column AR shows displacement of the 
final relaxed atomic positions from the refined diffraction- 
data sites, and the last column (xai) is occupancy by Al atom 
from Ref. [20] The experimentally determined structure has 
higher symmetry (space group Cmmm), and shows aver- 
aging over disorder manifested at M2 site, forming 0.8A- 
distant pairs of half-occupied positions. Occupancy of the 
M2 site couples with Al/Cu symmetry breaking of its adjacent 
M8a/b and M9a/b sites - all other sites marked a/b in the table 
(Sc3, M3, M10 and M13) have unique chemical assignement. 
Not surprisingly, largest displacement AR from the Cmmm 
diffraction-data model occurs for Al-occupied variants of the 
M8 and M9 sites (Cu is stronger X-ray scatterer). The av- 
erage Cmmm model can be therefore interpreted as thermal 
average over (locally) ordered patterns, in which flip of the 
M2 atom by ^0.8 A over pseudo-mirror plane perpendicular 
to c axis at z=Q is correlated with Al=Cu swap of adjacent 
M8a/M8b and M9a/M9b sites. 
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TABLE I: Fitted parameters for Al-Cu-Sc EOPP potentials. 
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4 


0.1821 


0.5 


-0.0072 


0.06 


0.78 


M12 


Al 


4 


0.5 


0.3050 


0.0034 


0.03 


0.69 


M13a 


Al 


8 


0.2201 


0.4187 


0.2728 


0.08 


0.93 


M13b 


Al 


8 


0.2194 


0.4199 


-0.2820 


0.07 


0.93 



this site has 0.5 occupancy in the refinement Ref.1201 



TABLE II: List of Wyckoff sites for oC104 structure minimizing 
EOPP energy, starting from random initial state, and using fixed ex- 
perimentally determined lattice parameters as the only input. Last 
two columns compare the refined structure with the Rietveld-refined 
structure of Ref.[20| see text. 



